Predicting Space Weather Effects on Close Approach Events 

Matthew D. Hejduk 

Astrorum Consulting 

Lauri K. Newman 

NASA Goddard Space Flight Center 

Rebecca L. Besser and Daniel A. Pachura 

Omitron, Inc. 


ABSTRACT 

The NASA Robotic Conjunction Assessment Risk Analysis (CARA) team sends ephemeris data to the Joint 
Space Operations Center (JSpOC) for conjunction assessment screening against the JSpOC high accuracy catalog 
and then assesses risk posed to protected assets from predicted close approaches. Since most spacecraft supported 
by the CARA team are located in LEO orbits, atmospheric drag is the primary source of state estimate uncertainty. 
Drag magnitude and uncertainty is directly governed by atmospheric density and thus space weather. At present the 
actual effect of space weather on atmospheric density cannot be accurately predicted because most atmospheric 
density models are empirical in nature, which do not perform well in prediction. The Jacchia-Bowman-HASDM 
2009 (JBH09) atmospheric density model used at the JSpOC employs a solar storm active compensation feature that 
predicts storm sizes and arrival times and thus the resulting neutral density alterations. With this feature, estimation 
errors can occur in either direction (i.e., over- or under-estimation of density and thus drag). Although the exact 
effect of a solar storm on atmospheric drag cannot be determined, one can explore the effects of JBH09 model error 
on conjuncting objects’ trajectories to determine if a conjunction is likely to become riskier, less risky, or pass 
unaffected. The CARA team has constructed a Space Weather Trade-Space tool that systematically alters the drag 
situation for the conjuncting objects and recalculates the probability of collision for each case to determine the range 
of possible effects on the collision risk. In addition to a review of the theory and the particulars of the tool, the 
different types of observed output will be explained, along with statistics of their frequency. 

1. INTRODUCTION 

The parameter at the core of the satellite Conjunction Assessment and Risk Analysis (CARA) problem is the 
probability of collision (Pc) parameter, which is the probability that two satellites will pass within a specified miss 
distance of each other. For most predicted close approaches, the dynamics of the situation allow the conjunction 
fundamentals to be projected into and analyzed in a plane perpendicular to the satellites’ relative velocity vector, 
permitting the calculation of the Pc using the “standard” formula given as the area integral shown below: 



in which r is the nominal miss distance miss vector between the satellites, C is the additive combination of the two 
objects’ covariance matrices, and A is the area representing the combined size of the two objects. [1] This Pc 
approach possesses a substantial advantage over the use of simple miss-distance calculations in that it incorporates 
the uncertainties of the two state estimates, so that for a fixed miss distance, as these estimates increase in 
uncertainty. uncertainties increase past a certain critical point the resultant Pc decreases and ultimately approaches 
zero. Of course, the rectitude of the calculation is governed by whether the objects’ covariances are complete 
(meaning that they attempt to take cognizance of all of the known error sources) and realistic (meaning that they 
reflect these errors properly and accurately). The question of covariance realism — whether object covariances are 
appropriately sized — is a general one that is now receiving broad attention within the conjunction analysis 
community, and a number of realism improvement initiatives are in various states of completion. [2,3] When all 
such efforts come to fruition, it is expected that the covariances supplied by the Joint Space Operations Center 
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(JSpOC), both in their circulated form and through realism enhancements that can be applied at the receiving end, 
will constitute a quite reasonable representation of the state estimate error. 

One area, however, that presently lacks robust error representation is active solar storm modeling. One reason 
this is true is that solar storm prediction modeling is a relatively new addition to routine space operations; as it 
became a routine feature of the JSpOC operational atmospheric density prediction model (JBH09) in the summer of 
2013, there has not been a large amount of subsequent historical solar storm data from which to evaluate the model’s 
behavior, and there have been no substantial studies of this to date. 1 Despite this lack of study emphasis, solar storm 
situations present particularly trying situations for CARA operations. Very large changes in atmospheric drag can 
be experienced between state estimate updates, and these changes can drive similarly large changes in the calculated 
Pc, stymieing the risk assessment process. In situations where some sort of active risk mitigation, such as a satellite 
maneuver, is being considered, the large changes in state estimates brought on by solar storms could actually make a 
post-remediation close approach situation more risky than if there were no remediation performed at all. In viewing 
the situation from the owner/operator’s point of view, what is needed is some sense of whether it is reasonable to 
take action (to remediate or elect not to remediate) based on the data at hand, or whether there is likely to be benefit 
to waiting for further updates, even though not acting right away will typically make the remediation process more 
rushed and frenetic. In short, some method is needed for determining how sensitive a particular conjunction is to 
changes in atmospheric drag so that the owner/operator can have a better sense of whether actions they are 
considering at any given moment are likely to be definitive or could become more dangerous. 

The present paper is organized in the following way. First, the general concepts of atmospheric drag and 
atmospheric density estimation will be treated, with special emphasis on the estimation problem in the presence of a 
solar storm. Next, the question of the sensitivity of conjunction assessment solutions to atmospheric density will be 
addressed and the Space Weather Trade-Space (SWTS) tool developed to explore this sensitivity will be discussed. 
Next, different “canonical" sensitivity response types will be described and shown, along with statistical data that 
document the relative frequencies of the different types. Finally, some concluding remarks will be offered on both 
the helpfulness and limitations of the tool, as well as suggestions for future work. 

2. SATELLITE DRAG ACCELERATION, ATMOSPHERIC DENSITY, AND SOLAR STORMS 

Satellite drag acceleration is the principal non-conservative force that affects low-Earth orbit (LEO) objects, 
which are defined as those satellites having a period less than 225 minutes. Drag acceleration is given by the 
relationship: 

1 C d A t/ 2 

a = — p— — V ( 2 ) 

2 M 

in which p is the atmospheric density, Cd is the drag coefficient, A is the satellite frontal area, M is the satellite mass, 
and V is the satellite velocity relative to the atmosphere. The group of terms CdAJM is called the ballistic coefficient 
and is often treated as a unit because for many objects, none of these terms is individually known; so the OD process 
treats them as a grouped unit and solves for this ballistic coefficient. To do this, as Eq. 2, makes clear, one must 
have estimated values for the satellite’s relative-to-atmosphere velocity and the expected density of that particular 
portion of the atmosphere. Since these three terms (atmospheric density, ballistic coefficient, and velocity) are all 
multiplied together, all of them can influence the atmospheric drag calculation substantially. It is also clear that, 
because they are part of a product, atmospheric density and ballistic coefficient can serve as proxies for each other: 
if one is curious about the effect on a satellite’s orbital dynamics if the atmospheric density were, say, twice as large, 
this could be investigated (in controlled circumstances) by making the ballistic coefficient twice as large. This 
aliasing between atmospheric density and the ballistic coefficient will become important in a later part of this paper 
in which dynamical sensitivity to atmospheric density will be explored. 


1 A broad investigation of the error behavior of the JBH09 model was recently executed and formed the basis for 
improvements to the JSpOC system’s consider parameter assignment to the drag variance term of its covariances. 
However, while this investigation did examine some solar storm data (i.e., situations when the disturbance storm 
time [Dst] index decreased below -75), the evaluation was a high-level treatment that did not consider how the errors 
evolved over the storm’s cycle and whether there was an estimation bias. 
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The process of estimating atmospheric density from empirical models, which is the form of model presently 
used in satellite operations, divides itself into two parts: one that models quiescent density (subject to routine solar, 
diurnal, and seasonal variations) and a second that models disturbed density, in which solar storm ejecta exert 
influence. In quiescent periods, nearly all of the atmosphere’s heating comes from absorption of extreme ultraviolet 
(EUV) radiation from the sun; so indices that reflect the amount of EUV received by the atmosphere can be fit 
through empirical models to actual density measurements and a reliable relationship established. To predict future 
atmospheric density values, one works from future predictions of the solar EUV indices; because these change 
slowly (with the eleven-year solar cycle), they can be and are predicted from autoregression fits of past data. 

Solar storms affect atmospheric density through a separate process called Joule heating; solar storm ejecta are 
mostly deflected away from the Earth by its magnetic field, but some of the material does enter through the polar 
regions (called the “polar cusp”) and directly interacts with the atmosphere, causing heating due to gas ionization 
and other molecular interactions. These heating activities are correlated with disturbances in the Earth’s magnetic 
field, measured by the older Kp/Ap geomagnetic indices and the newer disturbance storm time (Dst) index. The 
empirical modeling approach is similar to that for EUV heating in that empirical functional forms are used to 
establish a relationship between the observed geomagnetic indices and actual atmospheric density measurements, 
and the established correlation can be used to calculate density values using the geomagnetic indices as inputs. This 
approach, while functional, has always remained low-fidelity because the indices are always “chasing the action” of 
the solar storm and because the traditional models lack a real predictive value: the fact that a large solar storm has 
taken place can be known approximately ten minutes after eruption (the light-time from the sun to the Earth); but 
because it can take some 50 hours for the storm ejecta to reach the earth and begin to change the observed 
geomagnetic indices, satellite propagations more than 50 hours into the future will have extremely incorrect 
atmospheric density (and therefore drag) predictions past the point at which the storm ejecta reach the atmosphere. 

The JSpOC’s deployment of the JBH09 atmospheric density model has greatly ameliorated this situation by 
including a sub-model that attempts to predict solar storm timing and morphology from radiometric rather than 
geomagnetic indices, meaning that these predictions can be made 10 minutes (rather than 50 hours) after eruption 
and thus applied to satellite propagations essentially as soon as the storm has been detected. The particular sub- 
model incorporated into JBH09 is the Anemomilos model developed by Space Environment Technologies. [4] 

While the actual methodological operation of the model is complex, as an abbreviated summary it uses the observed 
(x-ray) magnitude of the flare to determine the mass of the ejecta and therefore the size and severity of the storm; it 
uses the intensity of the flare as a proxy for acceleration and thus the integral of this to determine flare velocity, 
which can be used to determine when the storm will arrive; and it uses the solar position of the storm (heliospheric 
latitude and longitude) to determine whether the storm will actually intersect the Earth and cause a disruption. From 
these values, the model selects a template that determines the Dst values versus time from storm initiation to 
conclusion, and these predicted values feed the Joule heating portion of the JBH09 modelling, which will thus 
generate atmospheric density estimates. 

The operational fielding of this sub-model is a major development in that routine space operations now possess 
a capability to predict and incorporate the expected effects of solar storms days before they actually reach Earth. 
Ironically, however, one is in a somewhat worse position in trying to determine the effect of the resultant change in 
atmosphere on the evolution of predicted satellite close approaches due to the potential for mismodeling of the solar 
storm. Without active storm compensation, one could be virtually certain that, during a solar storm, the actual 
amount of atmospheric drag experienced by the object would be greater than what the quiescent model would have 
predicted; one could then analyze the conjunction assuming an increased drag level to determine the probable effect 
on the dynamics. With active solar storm compensation automatically included in object trajectory predictions by 
the JSpOC when it propagates the object states and computes the time of closest approach, the overall magnitude of 
the drag modeling errors is typically expected to be smaller, 2 but one can no longer be sure of the manner of the 
mismodeling: the actual drag could now be either larger or smaller than the prediction, so the effect on the 
dynamics of both larger and smaller drag values must be assessed. This ambiguity of error direction complicates the 
error assessment process, and the lack of some sort of actual error function or model might lead one to conclude that 


2 But in a somewhat unexpected way this is not always true. In some cases the Anemomilos model will predict a 
storm that does not in fact materialize (it either turns out to be much less powerful than the predictions testified or 
the storm trajectory is incorrectly modeled and the storm ends up missing the Earth), and in such cases the model 
error is much larger than the predictions for the quiescent case. 
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really nothing of value can be done for model error assessment. On the contrary, given the manner in and timelines 
by which conjunction risk assessment decisions are made, there can be considerable value in determining the 
sensitivity of a particular event to atmospheric drag mismodeling. If the event is not sensitive to variations in the 
atmospheric drag, then the owner-operator can act with increased confidence on the information presently at hand, 
knowing that failures to model an upcoming solar storm correctly will not appreciably affect the conjunction risk 
assessment. If the event is in fact sensitive to the drag modeling, then the particular effects of mismodeling cannot 
be assessed; but if the former situation inheres in a considerable number of cases, then there is operational benefit to 
a tool that can assess this sensitivity. 

3. ATMOSPHERIC DRAG SENSITIVITY AND THE SPACE WEATHER TRADE SPACE (SWTS) 

As mentioned previously, Eq. 2 shows that the atmospheric density term and the ballistic coefficient are coupled 
as a product, meaning that the effect on a satellite’s dynamics of modifying one of these parameters can be emulated 
by a modulation of the other. It is relatively cumbersome to try to insert into a pre-packaged propagator an altered 
version of the output of an atmospheric density model, but it is quite straightforward to alter the satellite’s ballistic 
coefficient and, using this altered ballistic coefficient in propagation, examine the effect on the satellite’s position 
and the conjunction’s Pc. A simple analytical procedure thus suggests itself: to vary systematically the ballistic 
coefficients of the primary and secondary satellites of a conjunction, in a nested and gridded manner, to observe the 
effect on the overall Pc for each pair of modifications and to extract from this the sensitivity of the overall event to 
changes in the atmospheric density. It is important to vary the ballistic coefficients of the primary and secondary 
satellites independently because these two satellites may inhabit different orbit regimes and thus be differently 
affected by relative changes in the atmospheric density: a solar storm could have a large effect on one satellite’s 
dynamics and rather little effect on the other. 

While the concept itself is simple, there are subtleties of implementation that should be noted. First, given the 
more modest type of outcome desired by the analysis, namely the degree of sensitivity of the conjunction to 
atmospheric density mismodeling rather than any absolute evaluation of modeling error, it would at first seem to be 
acceptable to use any reasonable atmospheric density model for the investigation. However, it must be remembered 
that the JBH09 solar storm model is the only operational model presently able to make time-phased adjustments to 
the density to model a predicted storm that has yet to appear in the Ap/Kp indices; so it is only this same model that 
should be used to determine the sensitivity to atmospheric drag of an event that may occur during the solar storm 
compensation period — such an event may not even register at all with a model that does not contain this feature. 
Second, one must determine the degree to which the ballistic coefficient should be altered in both directions, and the 
granularity of alterations, in order to reveal the full morphology of the behavior and thus illuminate the trade-space. 
To get a sense of this, one should first explore the full range of atmospheric density values that can be produced by a 
robust atmospheric model. The output given in Figure 1 below was generated from the Jacchia 70 model for a 600 
km altitude at six different FlObar values (the 58-day moving average of the 10.7cm solar radio flux parameter) and 
the levels of F10 and Ap values shown on the x and y axes, respectively. Output total atmospheric densities are 
given by the contour colors as loglO density values (in kg/m 3 ). 

One can see that there is a total variation in density across all the response plots of a little more than two orders 
of magnitude, and even these probe the very extremes of the model’s response; so there would be no reason to vary 
the ballistic coefficient more than a total multiplicative range of an order of magnitude of increase and decrease. In 
experiments with different ranges of scaled ballistic coefficients, it was determined that a multiplicative variation 
range, in log space, of an increase and decrease of the ballistic coefficient of half an order of magnitude for both the 
primary and secondary was adequate to reveal the morphological features of the dependencies and allow 
categorization into a particular response type. 
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Fig. 1: Modeled atmospheric density values as a function of F10, FlObar, and Ap 


Therefore, to construct a tool that allows the determination of sensitivity of an event’s Pc to drag mismodelling, 
capitalizing on the coupling of atmospheric density and ballistic coefficient into one product, vary systematically the 
ballistic coefficients of the primary and secondary satellites of a conjunction within half of an order of magnitude, in 
a nested and gridded manner, to observe the effect on the overall Pc for each pair of modifications and to create a 
plot showing the sensitivity of the overall event to changes in the atmospheric density. CARA has followed this 
paradigm to build such a tool — the SWTS utility — to generate such sensitivity plots. 

4. SWTS RESPONSE CATEGORIES AND EXAMPLES 

The ability to generate a SWTS plot was implemented in CARA software and deployed operationally about 
eighteen months ago. A SWTS plot is produced whenever a conjunction event’s Pc exceeds IE-05; once this has 
happened, SWTS plots are generated for the event after each screening data update until the event time of closest 
approach has passed. Since this capability became operational, about 16,000 such plots have been produced. 
Examining this historical output, one observes that the results set can be categorized into three general types that 
further map into particular operational interpretations. A description of each of these categories, along with an 
example SWTS plot, is given in the following discussion. 

Fig. 2 illustrates an “on-ridge” SWTS plot, which is the most frequently encountered result category and, 
fortunately, the most operationally enlightening. In these plots, the nominal Pc is located at the center of the graph, 
and each axis gives the factor by which the primary or secondary satellite’s ballistic coefficient was shrunk or 
expanded. The contour colors represent the resultant Pc value for the encounter given the modified ballistic 
coefficients. In interpreting such a plot, what is important is not the actual calculated Pc values but the overall 
behavior of the contours and the location of the nominal Pc with respect to these countours. In an “on-ridge” 
situation, the nominal Pc is positioned on or very near (within half an order of magnitude in Pc space) the highest 
contour swath. This means that if the atmospheric density value predicted by the model is incorrect to some degree. 


5 


and thus the “true” Pc point is actually someplace on the graph other than in the very center, the Pc will be either the 
same as or lower than the nominal, but it will not be higher: the CA event risk will only decrease with atmospheric 
density model error. While of course it would be nice to have a model that lacked that error so that the 
owner/operator can make a remediation decision based on the true risk, in the above case he can act confidently on 
the information that he has, knowing that atmospheric density mismodeling will not make the situation worse than it 
appears presently. Of course, if new tracking data and a new OD set for primary and secondary are obtained, the 
entire solution may be different; and in such a case the Pc value may in fact increase. However, this increase is not 
driven by density mismodeling but by an improved knowledge of the actual situation fed by more recent tracking 
data. 
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Figs. 2 and 3: SWTS plot for an “on ridge” (above) and “flat” (below) situations 


Fig. 3 illustrates a “flat” situation, defined here as a plot with a total Pc variation less than an order of 
magnitude. Relatively uncommon, such a situation can be given essentially the same operational interpretation as an 
“on ridge” plot: atmospheric density mismodeling will not have a worrisome effect on the event risk; in fact, in this 
situation it will not have really any effect at all. In such a situation, one can counsel the owner/operator to make 
remediation decisions based on the information at hand, without additional worries about atmospheric density 
mismodeling. 

Fig. 4 gives an example of an “off-peak” situation, in which there is more than an order of magnitude’s total 
variation in the Pc throughout the plot and the nominal position is more than half an order of magnitude from the 
maximum value. Here, atmospheric density error could well increase the event risk, in some cases significantly. 
Density error could also decrease the risk. Because it is not known which effect would be present in any given 
situation, no useful recommendation can be given to an owner/operator in such a case; and thus here the SWTS is 
not helpful. However, if situations where it is helpful are prevalent, then the tool has operational value. 
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Fig. 4: SWTS plot for an “off-peak” situation 
5. SWTS -RELATED STATISTICS 

A question both of general interest and as an assessment of operational utility is the relative prevalence of the 
three different SWTS canonical types. The some 16,000 SWTS plots generated to date by the CARA operational 
tool were analyzed and categorized. CARA protected assets are grouped by orbit regime, and it makes sense to 
retain these same groupings (given in Table 1 below) for the present analysis. 

Table 1: CARA Orbit Regime Definitions 


Orbital 

Regime 

Definition 

LEO #1 

Perigee < 500 km & Eccentricity < 0.25 

LEO #2 

500 km < Perigee <750 km & Eccentricity < 0.25 

LEO #3 

750 km < Perigee < 1200 km & Eccentricity < 0.25 

LEO #4 

1200 km < Perigee < 2000 km & Eccentricity < 0.25 

HEO #1 

Perigee < 2000 km & Eccentricity > 0.25 
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MaximumPc/MinimumPc Order of Magnitude MaximumPc/NominalPc Order of Magnitude 


Fig. 5: CDF plots of MaxPc / MinPc and MaxPc / NominalPc 


Fig. 5 above shows two empirical cumulative distribution function (ECDF) plots of the SWTS database (with 
trivial results for which the maximum Pc = 0 removed from the sample) 3 showing the distribution of two different 
qualities: the ratio of maximum Pc to minimum Pc (left side) and maximum Pc to nominal Pc (right side). In each 
case, the base 10 logarithm of the ratios is what is plotted, so the x-axis renders the differences in “orders of 
magnitude.” In the left-hand graph, setting aside the results for LE04 and HEOl (which contain relatively small 
numbers of samples), one observes that about 80% of the cases have a dynamic range between maximum and 
minimum Pc of between four and seven orders of magnitude, meaning that most SWTS plots present a “ridged” 
response situation. It is also clear that only a very small percentage of the SWTS plots are of the “flat” designation: 
only a few percent show an overall dynamic range less than an order of magnitude. These statistics verify the 
operational intuition that most conjunctions are sensitive to atmospheric drag and therefore space weather effects; 
what remains to be determined is whether the typical case is a situation in which the nominal risk represents a level 
close to the maximum value. 

This question is addressed by the right-side graph in Fig. 5, which shows a CDF of the ratio between the 
maximum Pc in the SWTS plot and the nominal value, which appears at the center of the plot. There is more 
disparity among different orbit types here; but looking first at the response for regimes LEO 2 and LE03 (which 
generate the majority of the SWTS plots due to more protected assets residing in these regimes), sixty percent of the 
SWTS plots manifest less than half an order of magnitude’s difference between the maximum and nominal Pc 
values, indicating either a “flat” situation, which has already been shown to be rare, or a situation in which the 
nominal Pc is on a ridge. This latter case, shown here to be the majority result, is, as was stated previously, the most 
informative: here the owner/operator can be informed that the nominal risk represents an expected maximum; and if 
this maximum is below a level for which remediation is indicated, then the owner/operator can confidently embrace 
this conclusion despite any potential errors in solar storm compensation. 


3 As was mentioned earlier, SWTS plots are generated for any event whose Pc exceeds IE-05, and in such cases 
their generation is continued until the event passes the time of closest approach between the satellites; this practice is 
used because this event has been a significant event, and there is therefore interest in continuing to monitor it 
through its life cycle. In a significant number of cases, such events “fall off’ during the monitoring period, meaning 
that the risk drops to a level where the Pc is now smaller than the machine precision of the numerical integrator 
(around IE-324 for the MATLAB quadrature function) and thus set to 0. In these cases, the Pc is 0 for every grid 
point. In principle, these plots could be designated as “flat,” as they meet the strict definition for that category; but 
because they represent essentially discarded situations, it is a better practice to treat them similarly for these 
descriptive statistics. 
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These two results sets are nicely summarized in Fig. 6, which provides a stacked-bar chart of the relative 
composition, by response type, of each of the orbit regimes. The “Max Pc = 0” response is added here for context. 
One notes that, even with this category included, the “At Peak” situation is still a majority result for nearly all of the 
orbital regimes and, with the “Max Pc = 0” excluded, constitutes a supermajority. Given this frequency for the “At 
Peak” designation, one can ascribe a true utility to the SWTS: in the majority of cases, one can inform an 
owner/operator that the risk assessment information available represents a conservative evaluation and can inform a 
remediation decision regardless of any potential error in the JSpOC solar storm compensation model. 
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Fig. 6: SWTS category relative population by orbit regime 


5. CONCLUSIONS AND FUTURE WORK 

CARA has developed an operational software tool to examine the sensitivity of Pc to atmospheric density 
mismodeling, especially in the presence of detected solar storms. Results from this tool can inform owner/operators 
of situations in which the present Pc represents a maximum value with regard to atmospheric density modeling error 
and thus will not worsen due to density model updates. During the 18 months that this tool has been in use, statistics 
have been collected on observed categories of plot behavior. It has been observed that a large majority of the cases 
exhibit “on peak” behavior, which is the most helpful for operators to feel comfortable in making decisions about 
close approach mitigation in the presence of space weather events, thus confirming the tool’s utility. 
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